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Abstract. 

We discuss the possibility of a dimensional reduction of the Finstein equations in 
black-hole lattices. It was reported in previous literature that the evolution of spaces con¬ 
taining curves of local, discrete rotational and reflection Symmetry (LDRRS) can be carried 
out via a system of ODFs along these curves. However, 3-1-1 Numerical Relativity computa¬ 
tions demonstrate that this is not the case, and we show analytically that this is due to the 
presence of a tensorial quantity which is not suppressed by the symmetry. We calculate the 
term analytically, and verify numerically for an 8-black-hole lattice that it fully accounts for 
the anomalous results, and thus quantify its magnitude in this specific case. The presence of 
this term prevents the exact evolution of these spaces via previously-reported methods which 
do not involve a full 3-1-1 integration of Finstein’s equation. 
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1 black-hole lattices and dimensional reduction 

Black hole lattices, and in particular those conformally related to the 3-sphere have 
recently been the object of several studies [1-6] serving as toy models for the evolution 
of inhomogeneous universes [2, 3], for the propagation of gravitational waves in periodic 
spaces [4], for the exploration of cosmological models of non-trivial topology [5], and as an 
application of Regge calculus [6] . 

In [3], in particular, it was pointed out that the existence of a time-symmetric spatial 
hypersurface in these models, in addition to the high degree of spatial symmetry of the cell 
edges, implied that the evolution of the proper length of these subspaces was governed by 
a system of ordinary differential equations, and was thus decoupled from the surrounding 
spacetime. Here, however, we show that this is not the case: the symmetry is not sufficient 
to suppress one term, proportional to the curl of the magnetic part of the Weyl tensor, which 
contains spatial derivatives of the extrinsic curvature, and therefore prevents the reduction 
of the evolution equations to a simple, localized ODE system. 

In section 2, we present the arguments of [3] and the corresponding ODE system. In 
section 3, we integrate the Einstein equations numerically in 3-1-1 dimensions, and show that 
the comparison between the result and the solution of the ODE system features an anomaly 
which converges to a non-zero value in the continuum limit. In section 4, we rederive the 
equations of motion of black-hole lattices on the cell edges, and illustrate the origin of 
the curl term. We also show that in the 8-black-hole case, the additional term initially 
vanishes together with its first two time derivatives, but its third time derivative does not. 
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thus providing the final piece of evidence that the ODE system from [3] does not hold for the 
8 -black-hole lattice. Finally, we show in section 5 that the anomaly measured numerically in 
section 3 coincides with the curl term derived in section 4, and provide a fitting formula for 
representing this term analytically, so that the ODE system can still be used in combination 
with this source term. We provide some conclusions in section 6 . 

We will use the following index conventions throughout the paper: Greek indices /U, v, ..., 
will run from 0 to 3 and denote the spacetime objects. Objects defined on a spatial slice of 
dimension 3 will be denoted using Latin indices i, j, .... Tensorial quantities will be some¬ 
times written in the index notation {'Jij, Kij) and sometimes in the index-free notation using 
the bold typeface 7, K. 

2 Definition and properties of the LDRRS curves 

In this section we will introduce the mathematical formalism required to formulate the result 
of this paper. We will first restate the definition of a curve with a local, discrete rotational 
and reflection symmetry (LDRRS), discuss the effects of these symmetries and present the 
ODEs derived in [3]. The derivation of the reduced Einstein equations in [3] has been per¬ 
formed using the orthonormal frame approach in the version given by van Elst and Uggla 
[7]. This approach is less known than the standard ADM formalism and certainly less useful 
in numerical investigations where we have direct access to the 3-metric and the extrinsic 
curvature on a time slice, rather than the Ricci rotation coefficients or the commutation rela¬ 
tions of an orthonormal frame. We shall therefore present here the derivation of the reduced 
evolution equations from the ADM equations. Naturally, the final result does not depend on 
the formalism used. 

Consider the vacuum ADM equations in the normal gauge (corresponding to Gaussian 
normal coordinates, shift /3® = 0 , lapse a = 1 ) 

= -2Kij ( 2 . 1 ) 

Ay = Rij - 2Kik K'^j + K Ay (2.2) 

where 77 is the 3-metric on a spacelike hypersurface S, Ay denotes the extrinsic curvature 
of this hypersurface, A = A*^ its trace and Ay is the 3-dimensional Ricci tensor of 7 y. The 
constraint equations read 

R + K^ - Kij = 0 (2.3) 

(A*^’- A 7 *^).. = 0 (2.4) 

where R = R}^ and the covariant derivative is taken with respect to 7 y [8-10]. 

Following [3], we assume that on a given S there exists a curve A and a discrete group 
of symmetries G in the form of discrete n-fold rotations about A together with reflections 
through planes passing through A. More precisely, we assume that for each a £ G there exists 
a mapping Ra defined on a neighbourhood of A which preserves both the 3-metric and the 
extrinsic curvature: 


Kl = J ( 2 . 5 ) 

A:K = K, (2.6) 

A* denoting the pullback of a tensor by Ra- We assume it leaves every point in A invariant: 

Ra{p)=P iip£\. ( 2 . 7 ) 
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It follows from the assumptions above that Ra induces a mapping on the tangent space at 
every point p € \, i.e. -R* : TpS i—?■ TpS, which leaves both yjj and Kij at p invariant. 

We also assume that the action of G on TpS is the action of the group of discrete rotations 
and reflections. This assumption may be phrased in the following way: let r generate the 
rotations and let m and r generate the reflections. Then, in an appropriately chosen, properly 
oriented orthonormal frame Oj in TpS, we have 


i?:(ei) = ei 

(2.8) 

r,*/ \ . 27r 

Rr[^ 2 ) = COS - 02 - 1 - sm - 63 

n n 

(2.9) 

r,*/ \ • 2vr 27r 

Rri^s) = — sm - 02 - 1 - cos — 03 

n n 

(2.10) 

= ei 

(2.11) 

RUe 2 ) = -62 

(2.12) 

= es 

(2.13) 


where ei has been chosen to be tangent to A. Note that ei, 62 and 03 are not coordinate 
basis vectors. 

It is straightforward to see that the time development of A under the vacuum Einstein 
equations in normal coordinates will be a curve with local rotational and reflection symmetry. 
In [3] the authors prove that the assumption of invariance under (2.8)-(2.I3) restricts the 
form of vectors and tensors at points lying on A. In particular, a vector held X* invariant 
with respect to rotation (2.8)-(2.I0) has to be aligned along the curve A at every point p G X: 


Xp — Xi ei 


(2.14) 


and every rotation-invariant symmetric 2-tensor Sij = S'(jj) is a combination of the metric 
and a symmetric traceless tensor; 


5*. 

Sp = ^ 7 + 5'ii ( «! (8) «! - ^OL2 ® 0:2 - <8 0:3 


(2.15) 


where is the trace of S, Rn is the (1,1) component of Sij in the orthonormal frame e* 
and Q* is the dual co-frame of e*, i.e. a^{ej) = <58. On the other hand, every rotation- and 
rehection-invariant antisymmetric 2-tensor Aij = has to vanish at p: 

Ap = 0. (2.16) 


It follows quite easily from (2.14) that A must be a geodesic with respect to Indeed, 
let V be a tangent vector to A in any parametrization. The vector VvV is rotation-invariant 
since both the curve A and the metric are invariant as well. From (2.14) it must be 
proportional to ei, and thus also to v itself. After a suitable reparametrisation we obtain 

VvV = 0. 

In [3], following [7], a formalism was presented for simplifying the system (2.1)-(2.2) on 
a LDRRS curve. In this case, the only degrees of freedom of the metric tensor not suppressed 
by the symmetries are: 

ay = V7(Zi,Zi) (2.17) 

a_L = \/7 (Z2, Z2) = \/7 (Z3, Z3), (2.18) 
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where the vectors Zi, Z2 and Z3 are initially equal to ei, 62 and 63 respectively, and their 
components are assumed to be constant in time in the normal coordinate basis (Zi, Z2 and 
Z3 are the coordinate basis vectors if we choose the initial coordinate system appropriately). 
Note that the fnnctions ay and a±_ are snfhcient to reconstruct the metric tensor at any time: 

7 = (^1 (8) cui + 0^2 ® <^ 2 ) + a||(t)^cu3 (g) cj3, (2-19) 

where Ui is the dual co-frame of Z*, with constant components in a normal coordinate system 
basis. The functions evolve according to: 


ay 


( 2 . 20 ) 

ay 


a_L 

a_L 

= --E+. 

3 + 

( 2 . 21 ) 


E-y. is the only surviving component of the electric part of the Weyl tensor, given by 

E+ = -^Eije{e{ ( 2 . 22 ) 

= (2.23) 

being the normal to the constant time slice. According to [3] its evolution is likewise 
governed by an ODE: 

E+ = -3—E+ (2.24) 

a_L 

so that the evolution of the geometry on the curve is completely decoupled from its sur¬ 
roundings (in fact, the evolution of every single point on the curve is decoupled from all the 
others), and quantities that only depend on the metric tensor on the cnrve can be evolved 
using jnst the above system of ODEs (note that we will show in the following sections that 
(2.24) is missing an essential term which causes this decoupling to fail). 

Such a simplified scenario is particnlarly suitable for use as a numerical testbed, as one 
can compare the results of a full three-dimensional numerical evolution to the functions ay, 
a_L and Ey. defined above, and check to what extent the code reproduces the ODE system. 
We illustrate the result of this comparison in the next section. 

3 Numerical Relativity solution of an 5^ black-hole lattice spacetime 
3.1 Methods 

We solve the full 3-1-1 Einstein equations for an lattice nsing Numerical Relativity, 
allowing us to compute the metric everywhere, not just on the points of high symmetry. We 
use the open-source Einstein Toolkit [11] and Cactus [12] framework. We compute various 
lattice-related analysis qnantities using a Cactus code generated nsing Kranc [13, 14] and the 
xAct [15] tensor-manipulation package. Analysis of the nnmerical data was performed using 
SimulationTools for Mathematica [16]. 

We focus on the tesseract configuration, in which 8 identical black holes are arranged 
regnlarly on S^. To simplify the numerical treatment, we carry out the stereographic projec¬ 
tion, introduced in [ 2 ], from to B?, where one of the black holes, with bare mass mi = 4M 
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are 


is projected into the coordinate origin, another six, each with bare mass m2-7 = 4\/2A4, 
projected to x* = (±2,0,0)At, x* = (0,±2,0)At, x* = (0,0,±2)A^, and the eighth is pro¬ 
jected out to infinity (its presence, however, is revealed by an inner trapped surface at a 
coordinate radius of about 20A^). Here is a mass parameter, and the seven are 
the bare masses of the black holes in the Brill-Lindquist formalism (the asymmetry coming 
from the stereographic projection, as their physical masses are all equal to each other). The 
3-metric of the t = 0 surface is given by: 

~ '4^ ^ij 

where 5ij is the flat metric and V’ is the conformal factor, which takes the form: 

x,y,z) = 1 + 2^^ 

where rUn is the bare mass of the nth black hole and is the distance from it. The initial 
data is time-symmetric, so the extrinsic curvature is initially zero; 


(3.2) 


(3.1) 


Kij = 0 . 


(3.3) 


We will study one segment of the LDRRS curve corresponding to one of the edges of the 
lattice lying on the diagonal between the vertices x* = (2/3, 2/3, 2/3)W1 and x® = (2, 2, 2)W1. 
The angular coordinate (j) along this edge can be related to the Cartesian coordinates in 
the stereographic projection via 


<t> 


vr 

— — arccos 
4/1-1- sin 


3 V 1 — sin 


3^2 -4\ 
3^2 + 4 / 


(3.4) 

(3.5) 


where d is any of the coordinates x, y oi z. (j) = 0 {d = 2/\/3) corresponds to the midpoint 
of the edge, and most of our results will be presented using this point as an example. The 
vertices are located at </> = ±7r/6. 

We refer the reader to [2] for a discussion of this initial-data set, which also clarifies the 
role of the parameter A4 and the invariance of the spacetime under a rescaling of M. 

The spatial computational domain is |x®| < 24\/3A4 ~ 41.6A4, and we make use of the 
reflection symmetry in the x, y and z directions about the coordinate planes through the 
origin to restrict the domain to the octant where x > 0, y > 0 and z > 0. 

The Einstein equations are solved using the McLachlan [17] code with fourth order 
centred finite differencing for spatial derivatives and fourth order Runge-Kutta for the time 
integration. We integrate the Einstein equations in the BSSN formulation [18-20], a variant 
of the system (2.1)-(2.2). Eor more details, we refer the reader to [2]. 

In contrast to our previous work in [2], we do not use the standard binary black hole 
coordinate conditions (l-|-log slicing and gamma-driver shift), and then reslice the spacetime 
in postprocessing so as to use proper time as a time coordinate. Instead, we have recently 
found out that one can evolve this lattice directly in the normal gauge (unit lapse and zero 
shift), and still reach a proper time coordinate of t ~ 110A4 before the metric becomes 
degenerate in this coordinate system at the black holes, and further numerical evolution is 
not possible [21]. 
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In order to assess the effect of the numerical grid spacing on the results, we compared 
solutions with different overall grid spacings, labelled by the number of points, n, in one 
dimension in a certain coordinate distance. We report here the results for n = 32,40 and 
48. Since the space and time derivatives are computed with fourth order accuracy in the 
grid spacing h (x 1/n, and this is expected to be the dominant source of error, we expect the 
numerical error in the solution to scale as E = 0{h^) = 0(n“^). Hence, the errors at the 
three resolutions should be approximately in the ratio 1 : 0.41 : 0.20. 

Due to the different length scales in the system, we use mesh refinement to concentrate 
the computational grid points in regions where small length and time scales need to be re¬ 
solved, and avoid the prohibitive computational cost of using this same resolution everywhere. 
Mesh refinement is provided by the Carpet [22] code. The coarsest grid (level L = 0) has a 
grid spacing of ho = This was chosen so that the midpoint of the edge lies on a 

grid point at every resolution to avoid the need to interpolate data there. Refined regions 
are created around each black hole (BH) and around the edge where we wish to measure 
quantities accurately. Several levels of refinement are used, each level having half the grid 
spacing and time step of its parent region, resulting in a hierarchy of nested boxes with all 
regions on level L completely surrounded by regions of level L — 1. The BH at the origin, the 
six BHs at r = 2, and the edge are on levels L = 7, 6 and 5 respectively. The locations of the 
boundaries between refinement levels are found to have a significant effect on high-frequency 
numerical error measured on the edge. This was minimised by ensuring that the boundaries 
remain fixed in time and do not intersect the edge. 

3.2 Results and comparison with analytic ODE 

We now wish to determine whether the Numerical-Relativity (NR) spacetime satisfies the 
ODEs derived in [3] on the edge of the lattice. Using the NR spatial metric 7 ij, we compute 
a_L and ay as a function of time at the midpoint of the edge using (2.18) and (2.17). 

In the [3] solution, the evolution of a_L and ay is determined only by E_|_, and satisfies 
the system (2.20)-(2.21). By computing ay and d_L from the NR data, we determine E+ 
from both equations. The result is shown in figure la, and we see that the NR evolution is 
consistent with (2.18) and (2.17). 

According to [3], E_|_ also evolves according to the ODE given by (2.24) which involves 
a_i_ only, so that the system closes and the solution at the midpoint (as well as of any other 
point on the edges) decouples from its surroundings, as there are no spatial derivatives in the 
equation. In figure lb, we show the NR solution for E_|_ and —3^E_|_ which should agree if 
the ODE derived in [3] is correct. There is a significant disagreement. We have computed 
the numerical error bars for figure 1, but they are too small to be visible, indicating that 
the disagreement is a feature of the continuum Einstein equations, and not simply due to 
numerical error. Figure Ic shows a_L and ay computed from NR and compared with the 
solution from the ODE. There is a disagreement of up to 1% which is not accounted for by 
the relative numerical error of ~ 10“^. We thoroughly investigated all possible sources of 
error in the NR computation, and found none that could account for the discrepancy. 

For future reference, we define the anomaly A as the unknown additional term in the 
equation for E_|_. Hence, (2.24) from [3] gets modified to: 

E+ = -2>—E++A. (3.6) 

a_L 

In summary, the numerical results suggest that there is a term, A, missing in the evolution 
system of [3] which affects a± and ay at the level of 1% by t = 110A4. 
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tIM 

(a) Consistency of £’+ com¬ 
puted from the NR a_L and ay 



tIM 

(b) Violation of the ODE 
E+ = -3^E+ by the NR 
solution 



tIM 

(c) a_L and ay differ by up to 1% 
between the NR and ODE solu¬ 
tions. This is much larger than the 
relative numerical error of 10“^. 


Figure 1: Comparison between Numerical Relativity and ODE 


4 The evolution equations on a LDRRS curve 

We now turn our attention again to the ODE system formed by (2.20), (2.21) and (2.24), and 
in particular attempt to show that it is the reduction of the system (2.1)-(2.2) on a LDRRS 
curve. 

Equations (2.1)-(2.4) for -jij and Kij are not closed at a single point because of the 
presence of the spatial derivatives of yjj in the 3-dimensional Ricci tensor Rij. We will try 
to close the system by extending it to include Rij as a new evolution variable, for which we 
require the time derivative of Rij. Recall that the time derivative of the 3-dimensional Ricci 
tensor takes the form 


Rij — 2 {^^i\jk + j\ik - “ '^ij-,k^^ (^-1) 

where the indices have been raised by the inverse metric 7 *-^, i.e. 7 ^^ = 7 ^, 7 -^*'. We substitute 
( 2 . 1 ) to obtain 


4- = (4.2) 

We would like to get rid of the second derivatives in the equation above. This obviously 
requires permuting the indices in the expressions above. Swapping the first or the second 
pair of the indices is straightforward, but exchanging two indices between the two pairs is 
less obvious. We introduce the following notation for the antisymmetrisation in the indices 
2 and 3: 


Uijkl — Kij-,kl Rik-,jl- (^-3) 

After a tedious exercise in index manipulation we arrive at the following expression for the 
time derivative of the Ricci tensor 

% = (r4) 
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where || • || excludes an index from symmetrisation and R^j^i denotes the Riemann tensor of 
In three dimensions Rijki can be expressed entirely via the Ricci tensor and the Ricci 
scalar [23]: 

R-ijkl — Rjl '^ik Ril 'Ijk Rjk '^il Rik Tjfc 'Ijl 'lik^ ■ (^•^) 

If we substitnte the relation above into (4.4) we obtain 

Rij = - 2 + i?/ Kk^ + Kki (4.6) 

+\RKij + KRi, - ^RK^,j + [/(..)/, (4.7) 

i.e. the time derivative of the Ricci tensor expressed directly via 'jij, Kij and Rij and a single 
k 

additional term . The last term, i.e. the symmetrized contraction of Uijki, is now the 

only one involving the spatial derivatives of K^j and 'jij and thus not expressible directly via 
-jij, Kij and Rij. We will introduce a new notation for the symmetrised contraction: 

U„ = C,,/. (4.8) 

4.1 Properties of Uijki relation to the magnetic part of the Weyl tensor 

Before we proceed with the derivation of the reduced Einstein eqnations, we will discnss some 
of the properties of Uijki and elucidate its relation to the magnetic part of the Weyl tensor. 
Recall that the Weyl tensor in a 3+1 decomposition may be represented by its electric 

and magnetic parts defined via (2.23) and 

Bf,u = -^C'/xa/cA (4.9) 

respectively, being again the normal to the constant time slice, and denoting the 

totally antisymmetric volume form [24, 25] . Both tensors vanish in the normal direction and 
can be considered 3-dimensional, spatial objects. The magnetic part of the Weyl tensor in 
the ADM variables can be related to the tensorial curl of the extrinsic curvature: 

Bij = r]^^’‘Kikj. (4.10) 

We can easily prove that it is traceless and symmetric: first, we note that 

B\ = v^'^^Kik-i = 0 (4.11) 

becanse Kik is symmetric with respect to the exchange of the indices. The contraction of Bij 
with another volume form is equal to zero as well: 

B,j = - ( 7 ^* 7 'P - 7 ^P 7 '*) Kikj = - K’P = 0, (4.12) 

where the last expression vanishes becanse of the vector constraint equation (2.4). Note that 
in three dimensions this implies the vanishing of the whole antisymmetric part of Bij, so 
Bij = -B(ij)- 

The covariant derivative of Bij on the other hand can be related to U via 

Uijki = Bip.jrf (4.13) 
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The equation above, contracted with respect to the two last indices, yields 


Bip-i, (4.14) 

which has exactly the same structure as (4.10), i.e. Uij is proportional to the curl of Bij. 
Now, repeating the reasoning we have used for Bij above we may prove that the trace of 
(4.14) vanishes because of the symmetry of Bij} 

U\ = 0. (4.15) 

4.2 The reduced evolution equations 

Consider the tangent space at a point along a LDRRS curve A. We rewrite (2.1)-(2.2) and 
(4.7) assuming conditions (2.14)-(2.16) to hold and parametrizing the metric according to 
(2.17)-(2.18). We first note that the antisymmetric part of 17^.^ must vanish because of 
(2.16). Since it is also traceless it must be proportional to Un (see (2.15)). We obtain 
(2.20)-(2.21), where is the non-vanishing part of the electric Weyl tensor, but (2.24) now 
takes the form of 


=-3^i7+-|[/n, (4.16) 

with 17ii = U- e\e\ (notice that numeric indices always indicate frame components). In [3], 
the authors assume that this term vanishes due to the rotation and reflection invariance.^ 
We will show that this is not the case in general. As a result, we will identify this term with 
the anomaly A found numerically in section 3.2; 

Al = -^C/ii. (4.17) 

First let us consider the magnetic part of the Weyl tensor. Since Bij is composed of 
rotation-invariant rjijk and Kij-k it is rotation-invariant itself. Being additionally traceless 
and symmetric it must be proportional to Bn due to (2.15). From (4.10) we obtain 

Bn = A'i2;3 — A'i 3;2 = 2iFj[2;3] e) (4.18) 

Since ei is both rotation- and reflection-invariant, the last expression is the (2,3) component 
of a rotation- and reflection-invariant rank 2 antisymmetric tensor, so it must vanish at A 
because of (2.16). 

Now, since Uij is also traceless and is given by a very similar expression (4.14) to Bij, 
it would be tempting to repeat the argument above and conclude that l/n, together with 
the whole symmetric part of Uij, vanishes too. This would however be incorrect due to the 
following: unlike Kij appearing in (4.10), Bij in (4.14) is not reflection-invariant. Note that 
since its definition (4.10) involves the volume form r]ijk it changes its sign under reflections 
(2.11)-(2.13). Although Un can be put in a similar form to (4.18): 

Un = —2Ri[2;3] (4.19) 

^Note however that Uij does not have to be symmetric in general, unlike Bij. 

^It corresponds to the term proportional to e.y in equation (2.15) in the aforementioned paper. 

If it does not vanish then it appears later in the evolution equation (4.11). 
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we cannot now apply (2.16) because the antisymmetric 2-tensor in question e| is not 

reflection-invariant. The symmetry assumptions put no restrictions on the value of Un and 
there no reason whatsoever to assume that Uu vanishes identically along a LDRRS curve. 

We can give a simple and instructive analogy from the theory of electromagnetism and 
Maxwell’s equations. Consider a static configuration of electric and magnetic fields in a flat 
space exhibiting a similar rotation and reflection invariance with respect to a chosen axis. If 
the vector potential A is invariant then its curl needs to vanish at the symmetry axis and 
thus the magnetic field B = curl A = 0 along the axis due to the reflection symmetry, just 
like in the case of the magnetic Weyl tensor. But the curl of B does not need to vanish at the 
axis. Indeed, it is easy to create a configuration of the electromagnetic field in which there is 
a non-vanishing current flowing along the axis and thus j = curl B ^ 0. This is due to the 
fact that B, as a curl of a vector, is a pseudovector field, while j, which is a curl of a curl 
of a vector is again a regular vector field. The former must vanish because of the reflection 
symmetry, but the latter not. 

4.3 Uu and its time derivatives for the 8-black-hole initial data 

The initial data described in section 3.1 is time-symmetric, so the solution in normal coordi¬ 
nates satisfies 7 ij(f) = 7ij(—t)- R follows that the odd time derivatives of the metric ^^ 2 v +1 lij 
and of the Christoffel symbols ^ 2 n+i ^ = 0,1, 2 ..., vanish at f = 0 identically. So does 

the extrinsic curvature together with its even time derivatives From (4.3) and (4.8) 

we see that the same must hold for Uij, i.e. 

d2N 

Qf 2 N ~ 0 at f — 0 (4.20) 

and in particular Uij and Uij vanish initially. Direct computation reveals that for the initial 
data (3.3) and (3.1) the first derivative Uij = 0 along a LDRRS curve vanishes as well. 
The first non-vanishing time derivative turns out to be the third one. We have evaluated it 
as a combination of the partial derivatives of the conformal factor if. Since the expression 
involves up to 4th covariant derivatives of the Ricci tensor of the tensor manipulations 
and algebraic reduction were performed using Mathematica. The final result, considered 
along the LDRRS curve and after simplifications due to the symmetry, reads 

uil^ L=o = 2V’"'® (- 

+2 + 2V^(2,2,2) ^ ^(2,4,0)^ _ ^(6,0,0)^ ^5 

+ (^_8(V;(3A0))2 ^ ^(2,0,0) ('35^(4,0,0) _ 29 (^^(0,0,4) ^ ^^( 0 , 2 , 2 ) ^ ^(0,4,0)^^ 

^ g^(0,l,l) (^^(0,1,3) ^ ^(0,3,1)^ ^ 7^(1,0,1) ('^(1,0,3) ^ ^(1,2,1)^ 
+7^(14,0) (^^(1,1,2) ^ ^(1,3,0)^ ^ + 2 i/;(1’2,2) ^ ^(1,4,0) _ ^(5,0,0)^^^ ^4 

+ 18V,(0-04)^(0,1,0)^(0,1,3) ^ i3(^(0,0,1))2^(0,2,2) ^ ^3(^(0,1,0))2^(0,2,2) 

+ 18V,(0-04)^(0,1,0)^(0,3,1) ^ 2(V;(0-0,1))2^(0,4,0) ^ ii(i/>(0’1’0))2^(0,4,0)^ 
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(^7 (v^(0’0’3) ^ ^(0,2,1)^ ^(1,0,1) ^ 7 (^^(0,1,2) + ^(0,3,0)^ ^(1,1,0) 

+3^(0,0,1) ^^(1,0,3) +^(1,2,1)^ ^3^(0,1,0) (^^(1,1,2) ^^(1,3,0)^^ _g(^(0,0,l))2^(4,0,0) 

- 20 + 2i/.(0’2,2) ^ ^(0,4,0)^^^ ^3 

-6 (3 (V;(2.0,0))2 _ ^(1,0,0) 1^49 j'^(0,0,l) (^^(0,0,3) ^ ^(0,2,1)^ 

+^(0,1,0) (^^(0,1,2) ^ ^(0,3,0)^^ ^(1,0,0) _ 8 ('(^(0,0,1))2 ^ (^(0,1,0))2^ ^(3,0,0)^^ ^2 

-288 , (4.21) 

where we have assumed above that the first coordinate is aligned along the curve and x^, x^ 
are transversal. We have introduced here a short hand notation for the partial derivatives in 
the form of d{x^Y '^' Substituting the conformal factor 'll: from equation 

(3.2), we obtain that the numerical value of u[i\o) at the midpoint of the edge is 4.3 x 10“^^. 
4.4 Effect of Uii on the metric 

The addition of the term —3/2Uii to the ODE clearly affects the evolution of ay and a_L- 
We can estimate the effect by making a Taylor expansion of 0 ||(t) and a_L(t) about t = 0 and 
using the evolution equations (2.21), (2.20) and (4.16) to evaluate the Taylor coefficients at 
t = 0. We find that the effect of Uu appears first in the 0{t^) term. Using an overbar to 
represent the solution using the original ODE (2.24), i.e. without the Uu term, we find 

Ao|| = a, - a, =-^d?(0)«' + O(«*). (4.22) 

Aox = ax - ax = + 0(t*). (4.23) 

(4.24) 

where we have used the fact that a±, 0 || and E+ have the same value at t = 0 independent 
of the appearance of Uu in the ODE. 

At the midpoint of the edge, at t = 110A4, we find a relative error Aa||/a|| of about 
1% compatible with the NR results in figure Ic, and a relative error of 100% by t = 235A4. 
We conclude that the leading order contribution to Uu leads to a complete breakdown of 
the original ODE solution by this time, though we cannot determine whether higher order 
corrections are important here. 

5 Numerical Relativity calculation of Uu 

5.1 Consistency between NR and new analytical results 

In the previous section, we identified a term, —3/2Uu, ia the evolution equation for E+ which 
was assumed in [3] to vanish, but for which we find a nonvanishing third time derivative. We 
now aim to verify that the NR solution satisfies the new evolution equation (4.16), and that 
the numerically-non-zero anomaly A is indeed related to Uu by (4.17). Uu is computed in 
NR from covariant derivatives of the extrinsic curvature, 

Uii = -Kik./)e\e{. (5.1) 
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tIM 



tIM 


(a) Comparison of i?+with the original (b) Uu and —2/3A com- 


and new RHSs 


puted from NR 



tIM 

(c) A = 1 - {-2/3A)/Un 

computed from NR at sev¬ 
eral resolutions. 


Figure 2: NR solution demonstrating consistency with the new analytic results 


whereas A is defined via (3.6). 

Figure 2a shows a comparison between and the RHS of the original and new evolution 
equations. We see that the addition of the term —3/2Uii is necessary for agreement. In figure 
2b, we see that Un and —2/3A are found to be indistinguishable, and figure 2c shows that 
their relative difference, A = 1 — (—2/3^)/C/ii, converges to zero as the numerical resolution 
n is increased. The convergence is 4th order, as expected from the finite differencing order of 
the code. A exhibits high-frequency noise for t < 40AI which we attribute to error coming 
from the finite precision with which floating point numbers are represented in the code^. 
We have partially filtered the high frequency noise from the data in figure 2c to make the 
convergence more apparent. For t > 40A4, there are lower-frequency oscillations in the error 
which we attribute to numerical reflections from mesh refinement boundaries. 

For t > 20A4, at the highest resolution, we see that |A| < 3 x 10“^. Hence 

= 1.00000(3)1/11 (5.2) 

o 

in agreement with the analytic derivation in section 4.2. For t < 20A4, the ratio is still 
consistent with —2/3, but the relative error is larger since Un itself is small. 

We therefore see that the anomaly originally measured in the comparison of the 3-1-1 
Numerical Relativity results and the ODE system presented in [3] was due to the term Un 
derived in the previous section, but taken to vanish in the original derivation. 

5.2 Computation of Un and fitting formula 

We now present the NR computation of Un on the edge, and give a simple fitting formula 
for it that could be used along with (4.16) to solve the system via an ODE. 

Figure 3 is a contour plot of logio(t/iiA4^) as a function of t and (f), the proper time 
and the angular coordinate along the edge, respectively. The black solid and dashed lines 

^.4 depends on the third time derivative of a±, and an initial relative roundoff error of e ~ 10“^® with 
frequency uj ~ njAt ~ 80, for At the time spacing of output data points, will be amplified by a factor of 
a±/’a ±UJ^ when taking a third derivative, which leads to a relative error in 'a± comparable with that observed 
for the measured values of a± ~ 10^ and a ± ~ 10“^. 
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n=48 

n=40 

Fitting formula 


Figure 3: Contours of log^g Un computed using NR at two different resolutions (solid and 
dashed black contours) and from a simple polynomial formula fitted to the NR data (blue 
contours) as a function of t (time) and cj) (coordinate along the edge). The fit is performed 
in the region lAf < t < llOAf. 

represent contours of C/n computed at resolutions n = 48 and n = 40 respectively. For 
t > 10A4 the two resolutions are indistinguishable, indicating that the numerical error is 
small in comparison with Uu. For t < lOAf, there are regions, notably around (j) ~ —tt/S, 
where the numerical error dominates over Uu, which at these early times is 0(10“^®). The 
NR data satisfies 17ii(0, 4>) = ±7r/6) = 0 on the initial slice and the vertices as expected, 

since Uu = 0 there by symmetry. Note that the NR computation was performed in Cartesian 
coordinates {t,x,y,z) and has been transformed to (t, 0 = | — cos“^ ( 3 x^+ 4 )) plotting. 
While the spacetime is symmetric about </> = 0 in this is not the case in {t,x,y,z), 

hence the continuum solution is expected to show this symmetry in cj), but the numerical 
error is not. This is reflected in Figure 3. 

Since Uu is only available numerically, and appears to have a simple form in the regions 
in which it is well resolved, we provide a simple formula based on low-order polynomials in 
t and (p obtained via a least-squares fit to the NR data in the region 1A4 < t < 110A4, 
—7r/6 < (j) < 7r/6, corresponding to the edge of the lattice. The fitting formula is 

Uu = Y.Cp,{t/MYP'iM-^ P = 3,5,7 g = 0,2,4,6 (5.3) 

p,<i 

with the coefficients Cpq of {t/MY 0'^ given in Table 1. The error estimate in the last 
digit in parentheses is an indication of numerical truncation error. The contours of the fitting 
formula are shown in blue in Figure 3. 



6 8 12 24 24 12 8 6 
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cjY 


00 


{t/Mf 

7-040(1) 

xlO-i^ 

-3-558(6) X10-12 

4-23(4) xlO -12 

-2-30(7) xlO 

(t/Mf 

-3-810(2) 

xlO-i'^ 

2-361(9) xlO-i*^ 

-4-37(6) xlO-i® 

3-1(1) 

xlO 

{t/MY 

1-0965(9) X10-21 

-6-85(3) X10-21 

1-30(2) X10-20 

-9-8(3) 

xlO 


=12 

-16 

-21 


Table 1: Coefficients of {t/MY and in the fitting formula for Un determined from NR 


The region lAf < t < lOAi contains localised regions of high relative numerical error, 
but the small number of degrees of freedom in the fitting formula means that the fit is 
insensitive to these localised regions. The region t < in which the NR error dominates, 
is outside the fit region, and hence the fitting formula is an extrapolation in this region. 
For t > lOA^, 1^1 < vr/S, i.e. the regions where Ui\ is not close to zero, this fitting function 
approximates the NR result to within ±1%. In the regions t < lOA^ and |(?!)| > vr/8, the 
absolute agreement is within 

For t > lOA^, the NR and htting-formula curves are visually indistinguishable. 


5.3 Computation of 

We now wish to compute the third time derivative of C/n at </> = 0 from the NR data and 
compare with the analytic result obtained in (4.21). We cannot directly finite-difference the 
NR data near t = 0 because, as can be seen in figure 3, it is contaminated by numerical error. 
Instead, we compute the derivative by analytically differentiating the fitting formula. The 
fitting effectively averages out the very small numerical errors near t = 0 and uses information 
from f > 0 , where the errors are less significant, to obtain information about the derivative 
at f = 0 . 

The fitting formula (5.3) contains only a finite number of terms, so the coefficients 
cannot be directly identified with the coefficients in a Taylor series, and hence with the 
derivatives of Uu. However, as the number of terms in the fitting formula is increased, we 
expect the coefficients to approach the Taylor coefficients. We find that as both Pmax and 

Q'max are increased, C 30 appears to converge exponentially towards a limiting value. Taking 

(3) 

this to be the Taylor coefficient, we obtain an NR estimate for 17^' which can be directly 
compared with the analytic value obtained from (4.21): 


d^Un 

dt^ 


t=O,0=O 


4.3015(4) X 10 ^^j\4 ® Numerical 
4.30113 X 10-^^M-^ Analytic . 


(5.4) 


The NR error estimate in parentheses includes the effect of both numerical truncation error 
and of fitting using a finite number of terms, and we see that the NR derivative matches the 
analytical calculation within NR errors. We therefore have a high degree of confidence that 
the numerical solution and our understanding of the analytical system are correct. 


5.4 Effect of t /11 

In figure 4, we show the relative difference between ay computed from NR and from the 
original ODE, and compare it with the leading order analytic contribution computed in 
section 4.4 from the Taylor series. We see that the relative difference is dominated by the 
leading order term for as long as the NR computation lasts. We do not know whether this 
will continue past t = 110A4. 
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Figure 4: Comparison between the difference between the NR and ODE solutions for ay, 
and the leading order contribution computed from a Taylor expansion. 

6 Conclusions 

We have solved the full Einstein equations for an 8-black-hole lattice spacetime using 
Numerical-Relativity and found that the solution along certain LDRRS curves does not 
agree with the ODEs previously derived for this system. 

We therefore analysed the behaviour of LDRRS curves in vacuum spacetimes, and in 
particular showed that the evolution of certain symmetric subsets does not decouple from the 
surrounding spacetime, following a system of pure ODEs, as had been previously claimed. 
Instead, the variables ay, a± and E_|_, which capture all the metric and extrinsic-curvature 
degrees of freedom not suppressed by the symmetries, follow a system of ODEs with a source 
term Un, which itself depends on the spatial derivatives of the extrinsic curvature and can 
only, to our knowledge, be computed via Numerical Relativity. 

We have then computed this term using Numerical Relativity, both as an anomaly in 
the original ODE system, and via its expression in terms of the derivatives of Kij. These 
agree to within the relative numerical error of 3 x 10“®, strengthening our confidence in both 
the analytical study performed in sections 2 and 4 and the numerical infrastructure used in 
sections 3 and 5, as well as in [2]. 

Note that it is still possible to use the ODE system presented in [3], as long as one 
knows the source term independently. To this end, we provide a polynomial function, fitted 
from the NR data, in both proper time t and edge coordinate (j), which can be used up until 
t ~ llOAI. We have also computed the third time derivative of Un on the midpoint of the 
edge at t = 0 (which is the lowest non-zero time derivative of Uu initially), and compared it 
with an analytical derivation of the same quantity in the ADM formalism. These also agree 
within the relative numerical error of 10“^. 

We conclude by remarking that the evolution of the edges of black-hole lattices is not 
symmetric enough for a reduced-dimension calculation'^, and can, to our knowledge, only be 
performed via Numerical Relativity. In particular, the solution to the ODE system presented 
in [3] can be treated as an approximation, valid at early times, which may or may not be 
close to the true solution at late times. We have measured the error associated with such 
an approximation by direct comparison with a 3-1-1 integration of the Einstein equations, 

^Note that the vertices, on the other hand, possess enough symmetries for a full decoupling from the 
neighbouring points, leading to the solution = 0 and the 3-metric being constant. 
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and find it to grow to ~ 1% in the components of the spatial metric for t < llOA^. We 
have shown analytically that the leading order effect of the Un term is 0{t^) in the metric, 
and this is observed to a very good approximation in the NR results for t < llOA^. The 
duration of the NR computations presented here is limited by the use of the normal gauge, 
and at the present time, we have no way of assessing the error resulting from neglecting the 
Uii term at late times far from the time-symmetric hypersurface. We observe that if the 
growth were to continue, the metric would have 100% error by t ~ 235M. At late times, the 
results obtained from the system derived in [3] may well be qualitatively different to those 
that would be obtained by an evolution using the full, corrected, system. 
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